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I. INTRODUCTION 

Dilute gases of alkaline fermionic and bosonic atoms 
are superfiuid at very low temperature: Bose-Einstein 
condensates have been obtained in the case of bosonic 
atoms pj , while condensation of molecules (made out of 
two atoms) has been observed in the case of fermionic 
atoms ■ For fermionic atoms in the weakly interacting 
regime 1, where /cf is the Fermi momentum 

and a the s-wave scattering length) BCS superfluidity is 
expected in the case of attractive interatomic interaction 
(a < 0). 

A striking experimental evidence for BCS superfluidity 
is still missing, even though various signals which would 
be coherent with a superfiuid behavior have been ob- 
served in some experiments: the anisotropic expansion of 
the gas after releasing it from the trap |3|, the measure- 
ment of the gap |J] , the measurement of the frequencies 
and damping rates of the breathing modes [j| . 

However, the gap has been actually measured only in 
the strongly interacting regime and no experimental val- 
ues exist for the weakly interacting case. The anisotropic 
expansion on the one hand and the frequencies of the 
breathing modes on the other hand can be predicted 
within a hydrodynamic approach for a superfiuid gas 
[E @- In both cases the experimental observations 
agree very well with the hydrodynamic predictions, and 
this could actually be considered as an evidence for su- 
perfluidity. However, the predictions for a superfiuid gas 
are the same as those for a normal gas in the presence 
of collisions. It is true that at the very low temperatures 
achieved in these experiments Pauli principle is expected 
to inhibit collisions. However, the experimental measure- 
ments have been performed during the expansion of the 
gas after releasing it from the trap. In such a situation 
momentum space deformations are possible and collisions 
can survive even at very low temperatures. So far, it 
has not been possible to completely control this problem 
from an experimental point of view and, for this reason, 
no firm conclusions about superfluidity can actually be 
drawn. 

Another limitation is related to the hydrodynamic 
approach: hydrodynamics can be safely applied only 



within the limits of validity of semiclassical approaches, 
A 3> Ml, where A is the pairing gap and f2 is the trapping 
frequency. Effects from the finite size and inhomogeneity, 
governed by the finite trap frequency Vl, are neglected. 
Moreover, the hydrodynamic formalism has been devel- 
oped so far only for the case of zero temperature (T = 0) . 

In this article we deal with the excitation spectra in 
the normal and superfiuid phases of a dilute Fermi gas 
and we analyze how these spectra are affected by super- 
fluidity, both in hydrodynamic and microscopic descrip- 
tions. In order to study excitations similar to those ob- 
served experimentally (the breathing modes) we focused 
our attention on the monopole and quadrupole modes. 
However, while the breathing modes have been observed 
for a cigar-shaped gas (and the radial and axial frequen- 
cies have been measured), we restrict our analysis to a 
spherical gas for the sake of numerical tractability. More- 
over, while the experiment of Ref. [j| has been done for 
a strongly interacting gas, we treat a weakly interacting 
system. 

We analyze the excitation spectra within a finite- 
temperature mean-field approach which provides a mi- 
croscopic treatment for the system. The Bogoliubov- 
de Gennes (BdG) equations are solved for the 
ground state and the excitations are treated within the 
quasiparticle random-phase approximation (QRPA) |lfj| . 
This approach has already been developed for atomic 
Fermi gases in Ref. where the spin-dipole and the 
quadrupole modes have been analyzed. On the other 
hand, the monopole modes have already been studied 
and compared to a schematic model in Ref. 01 ■ 

In the present work we want to study systematically 
the effects related to the temperature and to the trap 
frequency of the system. In particular, we compare our 
results with the corresponding hydrodynamic ones in or- 
der to check the validity of the semiclassical approach. 
In addition to the strength distributions related to the 
excitation spectra, we also present the transition densi- 
ties which can give important information on nature of 
the collective modes. 

The article is organized as follows. In Sect. [H] we 
briefly sketch the quantum mechanical and semiclassical 
formalisms to describe collective modes in the superfiuid 
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phase and in the normal phase in the collisionless limit. 
In Sect. ITTTl results for the monopole and quadrupole ex- 
citations are shown: the dependence on the temperature 
and on the frequency of the trap are studied. In Sect. II VI 
we draw our conclusions. 



II. QUANTUM MECHANICAL AND 
SEMICLASSICAL FORMALISM 

In this section we will briefly review the theoretical de- 
scription of collective modes in trapped Fermi gases. As 
already mentioned in the introduction, one has to dis- 
tinguish between quantum mechanical ( "microscopic" ) 
and semiclassical approaches. The fully quantum me- 
chanical calculation consists in solving the QRPA equa- 
tions, which are the small-amplitude limit of the time- 
dependent BdG equations. At present such calculations 
are available only for systems containing up to ~ 10 4 
atoms in the case of a spherically symmetric trap. These 
conditions are quite far from the experimental ones, cor- 
responding to particle numbers of ~ 10 5 — 10 6 particles 
in a cigar-shaped trap. Up to now, the "realistic" condi- 
tions can only be treated within semiclassical approaches. 
The simplest semiclassical approach is the hydrodynamic 
theory. This theory is valid in the superfluid phase at 
zero temperature, since the pairing correlations keep the 
Fermi surface spherical during the collective motion of 
the system. However, hydrodynamics fails at non-zero 
temperature, unless the local equilibrium can be ensured 
by collisions. Since we are interested in the weakly inter- 
acting regime, the collision rate 1/r is very small com- 
pared to the frequency of the trap. In this "collisionless" 
regime, the Fermi surface becomes locally deformed dur- 
ing the collective oscillation. This cannot be described by 
hydrodynamics, but requires a description in the frame- 
work of the Vlasov equation. The latter is valid in the 
normal phase, i.e., above the critical temperature T c . In 
the intermediate temperature range < T < T c , a semi- 
classical theory is still missing. 



A. Quantum mechanical formalism (QRPA) 

The QRPA method has already been applied to 
trapped Fermi gases in the weakly [nj as well as in the 
strongly interacting regime |13| and here we will only give 
a short summary. 

We consider a gas of atoms with mass m in a spheri- 
cal harmonic trap with frequency f2, assuming that the 
atoms equally occupy two hyperfine states a = \,\. Be- 
cause of the low temperature and density of the gas, the 
interaction between the atoms can be chosen as a zero- 
range interaction and parametrized by the s-wave atom- 
atom scattering length a. In order to simplify the nota- 
tion, we will express all quantities in harmonic oscillator 
(HO) units, i.e., frequencies in units of Q, energies in 
units of Ml, temperatures in units of HQ/kg, and lengths 



in units of the oscillator length Iho — \Jfr/{rMl). Fur- 
thermore, instead of the scattering length we will use the 
coupling constant g = Ana jl ho as parameter of the in- 
teraction strength. 

As mentioned above, the QRPA describes small- 
amplitude oscillations around the equilibrium state 
within the BdG formalism. Therefore, the first step con- 
sists in solving the BdG equations [9j 

[Ho + W(r)]u n i m (r) + A(r)v n i m (r) = E ni 
A(r)u„ im (r) - [H + W(r)]v nlrn (r) = E nl v nim (v) 

for the static case. In this way we obtain a set of 
quasiparticle energies E n i and wave-functions u n i m and 
Vnim ■ In Eq. Q , Hq denotes the hamiltonian of the non- 
interacting HO minus the chemical potential [i, 

#o = i(-V 2 +r 2 )- M , (2) 

while the interaction is accounted for in a self-consistent 
way through the Hartree potential W and the pairing 
field A. Due to spherical symmetry, the wave functions 
can be written as 

u n im (r) = u n i {r)Y lm (0,(f>), (3) 
v n i m (r) = v n i{r)Yi m (9,(f>) . (4) 

The quantum numbers I and m are the angular momen- 
tum and its projection, while n numbers different states 
having the same I and m. In practice, the diagonaliza- 
tion of Eq. is done in a truncated harmonic oscillator 
basis, containing the eigenfunctions of the trapping po- 
tential up to a certain HO energy Ec = Nc + §, i.e., 

2{n - 1) + I < N c ■ (5) 

The self-consistency relates W and A to the wave func- 
tions u and v . The mean field W is just proportional to 
the density, i.e., 

W W = ffE ~ f( E m)) + u nl {r)f{E nl )} , 

nl 

(6) 

where 

m = Tzihr (7) 

denotes the Fermi function. The Hartree field is inde- 
pendent of the cutoff Nc if the latter is taken sufficiently 
large. The calculation of the pairing field A, however, 
is more complicated. The zero-range interaction leads to 
a divergence which in the case of uniform systems can 
be regularized in a standard way by renormalizing the 
scattering length. This regularization method has been 
generalized to the case of trapped systems by Bruun et 
al. 0] and developed further by Bulgac and Yu and 
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two of the authors As a result, the pairing field can 
be written as 



No 



A(r) = -g eff (r) £ ±?—u nl {r)v nl {r)[l - 2f{E nl )] 



ni 



4tt 



(8) 

with an effective coupling constant g e g which allows to 
include the contribution from states beyond the cutoff 
Nc within the Thomas-Fermi approximation (TFA). The 
explicit expression for g e g reads 

1 1 1 (k F {r) k c {r) + k F {r) \ 
hi- — — : — —~kc{r)J, 



9eff{r) g 2tt 2 



k c {r) - k F {r) 



(9) 

where k F and kc denote the local Fermi and cutoff mo- 
menta, respectively: 



k F {r) = a/2/x - r 2 - 2W(r) . 



fc c (r) = y/2N c + 3 - r 2 



(10) 
(11) 



Once the static BdG equations are solved, we can cal- 
culate the linear response of the system to a small time- 
dependent perturbation. Following Ref. E3> we have to 
compute the QRPA response function II, which is a 4 x 4 
matrix built out of 16 correlation functions: 



n(w,r/) = 



/«p T p T » «p T p;» «p T x» «PTX t ))\ 

«PiPT» «/S;Pi)) ((PIX)) ((PiX f )) 

«xp t » ((xPi)) «xx» ((xx f » 

\((X t /5T» ((X^l)) «*+*» ((x^)}/ 



(12) 



with the short-hand notation 
dt 



{{AB)) 



2tt 



([A(t,r),B(0,r')]), (13) 



where () means the thermal average. The operators of 
the normal and anomalous densities, p and x, are defined 
in terms of the field operators ip and ip* as follows: 



Pcr(t,r) = $l{t,r)i[)v{t,r) , 
X(t,r) = ^j.(t,r)Vi T (t,r) . 



(14) 
(15) 



In order to obtain II, we first compute the free or un- 
perturbed response function Ho, which is defined analo- 
gously to Eq. (|T^|> . but which does not include the effect 
of interactions between the quasiparticlcs. Thus IIo can 
be obtained by replacing the field operators ip in Eqs. 
JUJ) and lO by 

4>„(t,r) = ^2[b n i m<7 u n i m (r)e lEnlt 

nlm 

-^ nlm _y nlm (r)e- iE ^}, (16) 

where b and are annihilation and creation operators 
of non-interacting quasiparticles. Inserting the result- 
ing expressions into Eq. I|12|) and using the relations 



{b a , bp) = {bl bjj} = 0, {6 Q , b\) = S al3 (l - f(E a )), and 
{b\,bp) — f(E a )6 a /3, we obtain explicit expressions for 
the 16 functions contained in IIo in terms of the u and 
v functions and the quasiparticle energies obtained from 
Eq. JTJ. 

Due to the spherical symmetry of the trap and the 
rotational invariance of the interaction, excitations with 
different angular momenta do not mix. Therefore it is 
useful to decompose IIo into contributions of different 
angular momenta: 

no(«, r, r') = J2 n 0L (w, r, r')Y LM (0, <f>)Y^ M (9', 0') . 

LM 

(17) 

The QRPA response 11^ for angular momentum L can 
now be obtained from the quasiparticle response IIol by 
solving the Bethe-Salpeter integral equation 

II l (cj, r, r') = U ql (oj, r, r') 

poo 

+ / dr"r" 2 Jl 0L (u 1 ,r,r")GIl L (uj,r"y), (18) 
Jo 

where G accounts for the residual interaction between 
the quasiparticlcs: 



G 



/0 g 0\ 

g 

g 

\0 g 0/ 



(19) 



When calculating the 16 functions contained in Hol, 
one observes that two of them, namely those related to 
((x^x)) and ({jofi)), are divergent for Nc — > oo. This di- 
vergence has the same origin as that of the pairing field. 
Bruun and Mottelson E3 therefore suggested to use the 
same pseudopotential method as for the regularization of 
the pairing field in order to remove the divergence. How- 
ever, it is not clear how in their prescription, Eq. (7) in 
Ref. E3 > the contribution of states beyond the cutoff Nc 
can be approximated (as we did in the case of the pair- 
ing field by using the TFA), which is crucial for having 
convergence at reasonable values of the cutoff Nc- We 
therefore propose a simplified prescription: when calcu- 
lating n L, we have to restrict the sum to states below 
the cutoff, 2(ra — 1) + I < Nc- To compensate the re- 
sulting cutoff dependence, the interaction in the pairing 
channel must be replaced by the effective coupling con- 
stant given in Eq. ©. Thus, we replace G in Eq. 1181) by 
G e ff(r"), which is defined by 



G eff (r) = 



/0 g 

g 

g eff (r) 

\0 g eff (r) 



(20) 



One can show that, in the case of a uniform system, this 
simplified prescription coincides with the pseudopotential 
method in the limit of excitations with long wavelengths 
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and low frequencies. We have checked the convergence 
of the results using this regularization prescription. 

Finally, we have to say how physical quantities of inter- 
est can be extracted from the correlation function II. To 
that end it is useful to look at the spectral representation 

2>*M) = f*S SiLJ '^ [ , (2i) 

' — ' J U! — U! + l£ 

with 

S(o; ) r ) r') = --Vlm((p CT M) 

7T * — ' 

<J<j' 

= (1 - e-' T ) "—f-^ - E J + E i) 

ij 

xJ2MMr)\3)(j\P*'(r')\i), (22) 

aa' 

where \i) and \ j) are eigenstates of the many-body hamil- 
tonian with total energies Ei and Ej, respectively, and 
Z = J^i ex p{Ei/T). In the present QRPA formalism Eq. 
(1221 is evaluated using the four upper left elements of the 
II response function 1|12|) . obtained with Eq. (|18fl. 

In this paper we will consider excitation operators of 
the form 

yi(i,r)«r 2 y LJ tf(M)e-*' rt . (23) 

with L = (monopole excitations) and L = 2 
(quadrupole excitations). The corresponding strength 
function Sl(uj), which gives the excitation spectrum, is 
defined by 

/>OC poo 

S L (u)= drr A dr'r' 4 VSi(u;,r,r') ■ (24) 
Jo Jo 

Another interesting quantity is the transition density 
dp = p— po, where po denotes the density in equilibrium 
and p is the density of the excited system. In the case of 
zero temperature, where the stationary system is in the 
ground state |0), the transition density for u) = Ej — Eq 
is proportional to 

Sp(uj = Ej - Eo,t) ex £0'IMr)|0) • (25) 

cr 

In this case, the sum over i in Eq. (1221 reduces to one 
term (i = 0), and therefore the transition density can be 
obtained from 

[Sp(ui = Ej - E , r)] 2 oc / dw'S{uj', r, r) , (26) 

where S is supposed to be sufficiently small to avoid that 
other states than the selected one contribute. 



B. Superfluid hydrodynamics 

At zero temperature, superfluid hydrodynamics pro- 
vides the equations of motion for the density (per spin 
state) p(t, r) and the irrotational collective velocity field 
v(£, r) of the superfluid current (continuity and Euler 
equations) [l7|: 

p+V-(pv) = 0, (27) 

v = _ V f- + -^ + ^V (28) 
V 2 m m J 

These equations can equally be used for fermionic and 
bosonic systems, only the equation of state, relating the 
local chemical potential pi oc to the density p, must be 
adapted correspondingly. In the case of weakly inter- 
acting fcrmions, where the density can be regarded as 
independent of the pairing gap, this equation of state is 
given by the Thomas-Fermi relation 

, loc ( P ) = l^ +9P = L_g +gp . (29) 

In the static (equilibrium) case, Eq. (|28|l together 
with this equation of state gives immediately the usual 
Thomas- Fermi equation for the density profile po( r ) ; 

Woc[po(r)] + Vb(r) = /x, (30) 

which is valid in both the normal and the superfluid 
phase. While the TFA in the normal phase is valid if 
pioc is much larger than the discrete level spacing of the 
trapped system (HQ, in our case), superfluid hydrodynam- 
ics requires in addition that also the pairing gap A is large 
compared with the level spacing, which is much more dif- 
ficult to satisfy. 

Since the superfluid velocity field v is irrotational, it 
can be written as a gradient. In order to establish a 
connection with microscopic quantities, we write it in 
the form 

v(r) = -V^(r) . (31) 
m 

where if is related to the phase of the pairing field by 
A(r) = |A(r)|exp[2^(r)]. 

In this article we are interested in small-amplitude mo- 
tion. We therefore split the density and the external po- 
tential into their equilibrium values and small deviations, 
P = Po + dp and V ext = Vo + Vi, and expand Eqs. l(2T|) 
and l|28(l up to linear order in the deviations. In addition, 
as we did in the preceding subsection, we will specialize 
to the case of a spherically symmetric harmonic trap and 
use the corresponding HO units (h = m = f2 = 1), i.e., 
Vo = r 2 /2. We know that for an excitation of the type 
(|23|l the solution must be of the form 

<p(t, r) = <p(r)Y LM (0, 4>) expH^i) (32) 

and analogous for Sp. Furthermore, we are interested 
in the eigenmodes of the system, which persist even if 
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Vi = 0. Then Eqs. lj2*7j l and l|2*5|) can be transformed into 
an eigenvalue equation for the eigenfrequencies lo and the 
corresponding eigenfunctions ip(r), 



dpi, 



dp 



(^(rW)' -L(L + V 



-UJ 2 (P; 



(33) 



where /' means df /dr, and an equation for the transition 
density, 



, . /dfiioc 

Op = 1LJ 1 ; 

V dp 



V = 



(34) 



The numerical solution of Eq. I|33|l is not difficult. 
However, in the present article we are only interested 
in the lowest monopole (L = 0) and quadrupole (L = 2) 
modes. For these two modes, the velocity field v is practi- 
cally linear in r, and we can thus obtain a very accurate 
analytic approximation to the numerical solution. Let 
us start with the quadrupole mode (L = 2). We insert 



the ansatz ip 



into Eq. H33|l . multiply the equation 



by po(r) and integrate over d 3 r. By this integration the 
small deviations of the quadratic ansatz from the exact 
solution of Eq. l|3"3")l are averaged out and one thus ob- 
tains a very precise prediction for the frequency. After a 
lengthy calculation we reproduce the well-known result 



=2=V2, 



(35) 



which is independent of the interaction. 

In a similar way we can find an approximation for the 
eigenfrequency of the lowest monopole mode (L = 0). In 
this case the function tp has the form (p(r) w a — br 2 . 
Inserting this ansatz into Eq. i|33|) . taking the derivative 
with respect to r in order to get rid of the constant a, 
multiplying by r and proceeding in the same way as in 
the case of the quadrupole mode, we finally obtain 



ui L =o = 2i/l 



3Ei, 



8E, 



(36) 



pot 



where Ei nt and E pot are the interaction and potential 
energies, 



E int = J d 3 rgp 2 {v) 1 E pot = J d 3 rr 2 p {r). (37) 

Contrary to the quadrupole frequency, the monopole fre- 
quency depends on the interaction. Since Ei nt is nega- 
tive, the frequency lol=o is slightly lower than twice the 
trap frequency, 2f2. Finally, the ratio of the constants a 
and b, which is needed in order to compute the transition 
density Sp, can be determined from the condition that the 
integral over dp must vanish, since the total number of 
particles stays constant. 



C. Vlasov description 

Let us now consider a normal Fermi gas just above T c . 
In the weakly interacting limit, T c is very small as com- 
pared with the Fermi energy, i.e., except for the fact that 



the system is not superfluid, we can neglect temperature 
effects. We will also assume that the effect of collisions 
can be neglected. Under this condition the system cannot 
come to local equilibrium during the collective motion. 
In order to describe this effect, we will use the Wigner 
function f(t, r,p). In equilibrium and within the TFA, 
this function simply describes a Fermi sphere: 



/o(r,p) = e(p F (r)-p). 



(38) 



Out of equilibrium, if the particles do not undergo enough 
collisions to restore the isotropic momentum distribution, 
the local Fermi surface will assume a more complicated 
shape. The equation of motion for the Wigner function 
is the Vlasov equation 



/ = (W) • (V p /) 



(V r /), 



(39) 



where V(t,r) = V ex t(t,r) + gp(t,r) is the total 
(external+mean-field) potential and V r and V p are act- 
ing in coordinate and momentum space, respectively. 

Contrary to the hydrodynamic equations in the super- 
fluid phase, it is very difficult to solve the Vlasov equa- 
tion directly. We are therefore again looking for approx- 
imate solutions for the special case of small-amplitude 
monopole and quadrupole oscillations in a spherical har- 
monic trap. We will employ the "generalized scaling 
ansatz" [l8|, which has been used with great success 
to describe giant resonances in atomic nuclei and which 
has also been applied to trapped atomic Fermi gases [||. 
In this approach, the possible deformations of the local 
Fermi surface are restricted to quadrupolar shape. Intro- 
ducing a small displacement field £(t, r), one can write 



with 



/(t,r ! p) = / (r , ,p'), (40) 

r' = r-£(t,r), (41) 
p' = p-m£(i,r) + V r [p •£(*,!■)]. (42) 



The velocity field is then simply given by v = £, and the 
last term in Eq. I|42l) describes the deformation of the 
Fermi sphere. For the form of the velocity field we make 
the same ansatz as before, i.e., 



aVr 2 Y, 



LM\ 



(43) 



with L — (monopole mode) or L = 2 (quadrupole 
mode). In analogy to the procedure in the preceding 
subsection, we linearize the Vlasov equation l)39|) with 
respect to £, multiply by p • and integrate over d 3 p 
and d 3 r. Using Eqs. (|3TJ|l we reproduce after a tedious 
calculation the results originally derived in Ref. [||, 



WL=o 




3Eki 
4:E pot 



(44) 
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Note that the monopole mode has the same frequency 
in the normal phase as in the superfluid phase. This 
can be understood as follows. If the displacement field is 
purely radial (£ oc r), as it is the case for the monopole 
mode, one can see from Eq. I|40(l that the Fermi surface 
stays spherical. Therefore hydrodynamics gives the same 
frequency as the Vlasov equation. The frequency of the 
quadrupole mode in the normal phase, however, is higher 
than in the superfluid phase by a factor of approximately 
y/2. From Eq. gDJ one can see that in this case the Fermi 
surface gets a quadrupole deformation perpendicular to 
the deformation of the density profile in coordinate space. 
This deformation costs energy and therefore increases the 
frequency of the mode as compared to hydrodynamics. 



III. RESULTS 
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FIG. 1: Free quasiparticle response (dashed line) and QRPA 
response (solid line) of the monopole excitation as a function 
of the frequency u) (in units of the trap frequency f2) , for three 
different temperatures: ksT — 0, 1.4 hfl, and 3 HQ (from left 
to right). 



In this section we will compare QRPA and semiclas- 
sical results for monopole and quadrupole oscillations in 
a spherical trap. We are mainly interested in the limits 
of validity of superfluid hydrodynamics, since this theory 
is widely used in order to analyze experimental results. 
For instance, recent experiments showed that the axial 
breathing mode in a cigar-shaped trap follows the hydro- 
dynamic behavior throughout the BCS-BEC crossover, 
while the radial breathing mode deviates considerably 
from the hydrodynamic predictions |l9j |. especially on 
the BCS side of the crossover region. In this experi- 
ment the system was still very strongly interacting even 
on the BCS side of the crossover (the strongest devia- 
tions happened when fcp-|a| was of the order of 2), such 
that our weak-coupling theory (valid for kp\a\ <C 1) can- 
not directly be compared to that experiment. Neverthe- 
less, it is clear that the limits of validity of hydrodynam- 
ics should be clarified. It is known that hydrodynamics 
works at zero temperature and if the level spacing Ml is 
much smaller than the gap A, but both conditions are 
generally not fulfilled in the experiments. Since exper- 
iments cannot be done at zero temperature, it is inter- 
esting to see what kind of temperature effects can arise 
below the critical temperature T c . The second condition 
is also very strong, especially if the trap is strongly de- 
formed and the transverse trap frequency is large, and it 
is therefore important to know up to which ratio Ml/ A 
hydrodynamics can be trusted. 
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FIG. 2: Free quasiparticle response (dashed line) and QRPA 
response (solid line) of the quadrupole excitation as a function 
of the frequency ui (in units of the trap frequency Q.) for three 
different temperatures: fcsT = 0, 1.4ftf2, and 3h£l (from left 
to right). 



zero temperature hydrodynamics should work very well. 

In Figs.Q]and[21we show the monopole and quadrupole 
response functions, respectively, for three different values 
of the temperature. The figures on the left show the re- 
sponse at zero temperature. The solid lines correspond 
to the QRPA results while the dashed lines represent the 
free quasiparticle response. In principle, the response 
function consists of a very large number of discrete lev- 
els. For the purpose of graphical presentation, these delta 
functions must be smeared out, and we therefore intro- 



A. Temperature dependence 

In this subsection we will study how the properties of 
collective modes change in the small temperature range 
from zero to the critical temperature T c . For this inves- 
tigation we are using the parameter set fi — 32 Ml and 
g = —0.965 (in HO units). With these parameters, the 
number of particles is approximately 17000 and the gap 
in the center of the trap at zero temperature is approx- 
imately 6 M7; one can therefore expect that at least at 





T -- 
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T>T C 




QRPA 


hydro. 


(Q)RPA 


Vlasov 


L = 


1.9 


1.88 


1.9 


1.88 


L — 2 


1.4 


V2 


2.2 


2.22 



TABLE I: Frequencies (in units of the trap frequency fi) of 
monopole (L = 0) and quadrupole (L = 2) modes for /i = 32 
hfl and g — —0.965 (in HO units) at zero temperature and 
above T c . The QRPA results for T > T c were obtained with 
T = 3m/k B - 
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FIG. 3: Transition densities for the collective monopole (left 
panel) and quadrupole (right panel) modes as a function of 
r (in units of the oscillator length Iho), at T = 0. Solid 
and dashed lines represent the QRPA and the semiclassical 
results, respectively. 



duce a small imaginary part of e = 0.015 il in the de- 
nominators of the correlation functions [see Eq. I|22l) ]. 
For T=0, the QRPA quadrupole response shows one sin- 
gle collective peak whose frequency is very close to that 
predicted by hydrodynamics (see Table I). The QRPA 
response is completely different from the free quasiparti- 
cle response, which has a broad and almost continuous 
distribution of strength between ~ 1.8 Q and ~ 2.7 fl. 
As has been realized before 0, the threshold of 
the two-quasiparticle strength is related to the energy 
of the lowest-lying quasiparticles which are located near 
the surface of the atomic cloud. 

In the case of the monopole mode the good agree- 
ment between QRPA and hydrodynamics (Table I) is 
even more surprising than in the case of the quadrupole 
mode, since the frequency of the monopole mode is so 
high that it lies in the two-quasiparticle continuum (see 
dashed line in Fig. ^| and one would therefore expect a 
certain amount of Landau damping. 

Apart from the study of the frequencies of the col- 
lective modes, the comparison between hydrodynamics 
and QRPA can be extended also to the analysis of the 
character of such modes. We display in Fig. the tran- 
sition densities of the two collective modes, which, since 
the density profile is known, can be related to the veloc- 
ity field [see Eq. pij)]. The normalization of the QRPA 
transition density is obtained from the integral of the cor- 
responding peak in the strength function, while that of 
the semiclassical transition density has been adjusted to 
the QRPA one. We see that the simple formulas from 
Sect. Ill 51 are in good agreement with the QRPA transi- 
tion densities. However, the QRPA transition densities 
exhibit small Friedel-like oscillations, especially near the 
surface where the gap is small and the local Fermi surface 
is therefore relatively sharp. 

Let us now consider an intermediate temperature be- 
tween and T c . For the present set of parameters 
the critical temperature is T c w 2.8 hfl/ks', we there- 
fore choose T = lAMl/ks ~ T c /2. As can be seen 



in the middle of Figs. ^ and |2 due to the presence of 
thermally excited quasiparticles the free quasiparticle re- 
sponse starts now already at u> — 0. As a consequence, 
both the collective monopole and quadrupole modes be- 
come strongly fragmented and damped. Qualitatively, 
this strong Landau damping at temperatures between 
zero and T c could be related to the damping mecha- 
nism which is responsible for the experimentally observed 
damping of breathing modes on the BCS side of the BEC- 
BCS crossover 19]. Interesting is also the double-peak 
structure which can be seen in the quadrupole response, 
as if there were two damped modes, one corresponding to 
the hydrodynamic mode and another one corresponding 
to the quadrupole mode in the collisionless normal phase 
(see below). This can be interpreted in the sense of the 
two- fluid model |20j , which states that between T = 
and T = T c the system effectively behaves as if it con- 
sisted of a mixture of normal and superfluid components. 

Now we increase the temperature further to T = 
ZMl/ks, which lies slightly above T c , i.e., the system 
reaches the normal phase, but still the temperature is 
very low compared with the Fermi energy. In the normal 
phase, the BdG equations become identical to the usual 
Hartree-Fock equations, and the QRPA becomes equal to 
the usual RPA. In the case of the monopole mode (right 
panel of Fig. QJ, the QRPA response is almost identical to 
that at zero temperature (left panel of Fig. ^) , although 
the free quasiparticle response is quite different. Again 
there is one collective mode having the same frequency 
as at T = 0. This is not very surprising. As mentioned 
in the preceding section, the Vlasov equation predicts 
the same frequency as superfluid hydrodynamics, since 
in the case of the monopole mode there is no deforma- 
tion of the local Fermi surface. This is different in the 
case of the quadrupole mode (right panel of Fig.0). Also 
here a collective mode reappears, but it is situated at a 
different frequency than at zero temperature. The higher 
frequency in the normal phase compared with the super- 
fluid phase is due to the Fermi surface deformation and 
is well described by the Vlasov equation (cf. Tabic QJ. 



B. Dependence on the size of the system 

Let us now investigate the importance of the discrete 
level spacing at zero temperature. In the case without 
superfluidity, the semiclassical h — * limit (TEA in equi- 
librium and the Vlasov equation in the dynamical case) is 
known to work very well if the number of particles is suffi- 
ciently large. The reason is very simple: The only dimen- 
sionless parameter on which corrections can depend is 
Ml/ (a, which becomes very small for large numbers of par- 
ticles. In the current experiments involving ~ 10 5 — 10 6 
atoms this type of corrections is completely negligible. 
For our study we choose, as in the preceding subsection, 
a chemical potential of \i — 32 Ml. This is large enough 
to make these corrections small, and the numerical cal- 
culations are still tractable. The corresponding numbers 
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a 


N 


A(0) 


-0.965 


16500 


6.0 


-0.8 


15000 


2.9 


-0.7 


14300 


1.4 


-0.636 


13900 


0.7 



TABLE II: Chosen values of the coupling constant g (first 
column; in HO units) and corresponding results for the num- 
ber of particles, N (second column), and for the gap at the 
center of the trap, A(0) (third column; in units of HQ). The 
remaining parameters were fixed to pi = 32 hfl and T = 0. 
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FIG. 4: Unperturbed response (dashed line) and QRPA re- 
sponse (solid line) of the quadrupole excitation as a function 
of the frequency ui (in units of the trap frequency £1) for T = 
and \x = 32 hfl and four different values of the coupling con- 
stant: g = -0.965, g = -0.8, g = -0.7, and g = -0.636 (in 
HO units; from top to bottom). 



of atoms lie between ~ 14000 and ~ 17000 depending on 
the chosen values of the coupling constant g due to the 
Hartree field (see Table HI)). 

In the case of superfluidity, however, another dimen- 
sionless parameter becomes important, which is Ml/ A. 
Since in the BCS phase A <C /U, this parameter is not 
necessarily small even if the number of particles is very 
large. In order to study the validity of hydrodynamics as 
a function of Ml/ A, we change A by varying the coupling 
constant g between —0.636 and —0.965 (in HO units). As 
a measure for A we take its value at the center of the trap, 
A(0). The values of A(0) corresponding to the different 
coupling constants are listed in Table ITT1 

We are now going to analyze the finite-size effects on 
the quadrupole response function by using the different 
values of the coupling constant listed in Table |HJ Note 
that, since we are using HO units, changing the cou- 
pling constant g cx a/ 1 ho is equivalent to changing the 
oscillator length Iho and thus the radius of the cloud 
R = \J1\i/Ml Iho- Anyway, as argued above, the impor- 
tant parameter for finite-size effects is the ratio Ml/A(0) 
and not the cloud size itself. 

For the strongest coupling, g = —0.965 (in HO units), 
the central value of the gap, A(0), is large compared with 
Ml, and hydrodynamics works almost perfectly at zero 
temperature, as we have already seen in the preceding 



subsection. Fig. 0] shows, from top to bottom, the evo- 
lution of the quadrupole response at T = for decreas- 
ing coupling constant g, i.e., for increasing importance of 
the discrete level spacing. Besides the QRPA response 
(solid lines), we also show the free quasiparticle response 
(dashed lines). For g = —0.8 (in HO units), the gap at 
the center is still larger than Ml by a factor of three, 
but now we find considerable deviations of the QRPA 
response from the hydrodynamic result. Since the free 
quasiparticle response is now shifted to lower frequencies, 
the hydrodynamic mode becomes fragmented, which ex- 
perimentally would show up as damping effect, and its 
frequency (u> ~ 1.1 fi) lies below the hydrodynamic pre- 
diction (V2fi). For g = -0.7 and g = -0.636 (in HO 
units), the central value of the gap is comparable to Ml 
and it is clear that hydrodynamics must fail. Indeed, 
the QRPA response becomes more and more similar to 
the free quasiparticle response which in the case of weak 
pairing looks very different from the strong-pairing case. 
The double-peak structure is a consequence of the two 
types of transitions which are allowed by the selection 
rules of the harmonic oscillator, i.e., transitions inside an 
oscillator shell (SN = 0, where N denotes the number of 
oscillator quanta) and transitions with 5N — 2. As the 
interaction decreases, the strength of the 6N — tran- 
sitions becomes less important while the 5N = 2 transi- 
tions become stronger. This can be understood from the 
fact that in the limit of a noninteracting harmonic oscil- 
lator without pairing (g — ► 0) the SN = transitions are 
forbidden by Pauli principle and only the SN = 2 transi- 
tions survive. In this limit the response has a single peak 
at oj = 2f2, in exact agreement with the prediction from 
the Vlasov equation |23j. In the semiclassical language, 
one can say that in this case the pairing is too weak to 
restore the spherical shape of the Fermi sphere during 
the oscillation, and therefore one finds the normal colli- 
sionless frequency instead of the hydrodynamical one. 



IV. SUMMARY AND CONCLUSIONS 

In this article we have studied the properties of collec- 
tive monopole and quadrupole modes in superfluid Fermi 
gases in the BCS phase (fcp|a| <C 1, a < 0) in a spheri- 
cal harmonic trap. Having briefly recalled the quantum 
mechanical and semiclassical formalisms (QRPA, hydro- 
dynamics, Vlasov equation), we presented numerical re- 
sults and compared the different formalisms. Our main 
interest was focused on two types of effects: temperature 
and finite-size effects. Both cannot be treated within 
the semiclassical approaches available in the present lit- 
erature, and they can therefore only be studied in the 
framework of the fully microscopic QRPA formalism. 

In the case of a sufficiently large system (large meaning 
A ^> Ml), superfluid hydrodynamics can be used to de- 
scribe the properties of collective modes at zero tempera- 
ture. Our results confirm earlier findings |ll| which show 
that already for parameters which lead to A(0) = &MI 
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the extremely simple theory of superfluid hydrodynam- 
ics is in almost perfect agreement with the numerically 
heavy QRPA method. This is not only true for the fre- 
quencies, but also for the transition densities, i.e., the 
velocity fields associated with the collective modes. How- 
ever, experiments can never be done at zero temperature. 
The critical temperature T c being extremely low, it is 
clear that already at very low temperatures between 
and T c the properties of the collective modes must un- 
dergo dramatic changes. This is evident if the hydrody- 
namic frequency (T = 0) is different from that in the 
collisionless normal phase (T = T c ), like in the case of 
the quadrupole mode. In the case of the monopole mode 
we also find a strong temperature dependence, although 
its frequency at T — is the same as at T — T c . In the 
intermediate temperature range between and T c the col- 
lective modes exhibit strong Landau damping. When the 
critical temperature is reached, the damping disappears 
and the collective modes can be very well described by 
the semiclassical Vlasov equation within the generalized 
scaling approximation. 

It is interesting to compare these temperature effects 
with those found previously in the case of the twist mode 
|2l| . which is an excitation where the upper hemisphere 
rotates against the lower one. Near T c , the behavior is 
rather similar: At T = T c the twist mode is a collective 
mode which can be described by the generalized scal- 
ing approximation to the Vlasov equation and whose fre- 
quency is slightly higher than the trap frequency. If the 
temperature is lowered, the twist mode becomes strongly 
damped, like the quadrupole and monopole modes. How- 
ever, an important qualitative difference appears near 
zero temperature. Since the velocity field of the twist 
mode cannot be written as a gradient, the twist mode 



disappears completely at zero temperature, whereas the 
quadrupole and monopole modes have an irrotational ve- 
locity field and they reappear at zero temperature as hy- 
drodynamic modes. In the case of the twist mode, the 
disappearance of the 1/to weighted integrated strength 
could be well described within a rather simple two-fluid 
model 0, |22| • It remains to be studied if a generaliza- 
tion of the two-fluid model to the dynamical case can also 
explain the damping of the quadrupole and monopole 
modes and the two-peak structure in the quadrupole re- 
sponse function at temperatures between and T c . 

In addition to temperature effects, we studied how the 
properties of the quadrupole mode change at zero tem- 
perature when the condition for the validity of the hy- 
drodynamic approach, A ^> Ml, is no longer satisfied. 
For parameters leading to A(0) « 'iMl the QRPA al- 
ready shows considerable deviations from the hydrody- 
namic theory. In the case of the quadrupole mode, the 
frequency for these parameters is found to be lower by 
20% than the hydrodynamic prediction, and a certain 
fragmentation of the excitation spectrum (i.e., damping 
of the collective mode) can be observed. If A(0) « Ml, 
the hydrodynamic mode has more or less disappeared. 
At the same time, a fragmented strength appears in the 
excitation spectrum near the frequency of the collective 
quadrupole mode in the normal collisionless phase. 

These results should be kept in mind when frequen- 
cies of collective modes measured in experiments with 
strongly deformed traps are compared with the hydrody- 
namic predictions. Due to the strong deformation, the 
radial trap frequency fl r is often much higher than the 
axial one, f2 z . Even in the case of strong pairing, the gap 
might be of the order of, say, 2>Ml z , and considerable 
deviations from hydrodynamics are possible. 
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